package model.util;

public class QuickFourierTransform {
	private static final int FFTDHMAXSIZE = 32768;

	public static void fftdhts(Compl[] x) {

		/*int len = x.length;
		int step = 0;
		if (len % 2 == 0) {
			while (len != 1) {
				if (len % 2 != 0) {
					throw new IllegalArgumentException(
							"Размер области не является степенью двойки");
				}
				len /= 2;
				step++;

			}
		}*/
		int i,j,S=x.length/2;

		Compl[] tab = new Compl[x.length];

		for(j=0;j<=S;j++)

		  {tab[j] = new Compl();tab[j].re=Math.cos(Math.PI*j/S); tab[j].im=-1*Math.sin(Math.PI*j/S);}
		fftdhts(x, x.length, tab);
	}

	private static void fftdhts(Compl[] x, int Size, Compl[] tab) {

		if (Size > FFTDHMAXSIZE || Size <= 0) {

			throw new IllegalArgumentException(
					"Size must be <= FFTDHMAXSIZE and >= 0");
		}

		int k = 1, Pow2 = 0;

		while (k < Size) {
			k *= 2;
			Pow2++;
		}

		if (k != Size)
			throw new IllegalArgumentException("Size must be power of 2");

		for (int i = 1; i < Size - 1; i += 2) {
			x[i].re = -x[i].re;
			x[i].im = -x[i].im;
		}

		int j = 0;

		for (int i = 0; i < Size - 1; i++) {

			Compl t;

			if (i < j) {
				t = x[j];
				x[j] = x[i];
				x[i] = t;
			}

			k = Size / 2;

			while (k <= j) {
				j -= k;
				k /= 2;
			}
			j += k;
		}

		for (int l = 0, ll = 1; l < Pow2; l++, ll *= 2) {
			for (j = 0; j < ll; j++) {

				Compl tv = new Compl();
				tv = tab[j << (Pow2 - 1 - l)];

				for (int i = j; i < Size; i += 2 * ll) {
					Compl t = new Compl();
					k = i + ll;

					t.re = x[k].re * tv.re - x[k].im * tv.im;
					t.im = x[k].re * tv.im + x[k].im * tv.re;

					x[k].re = x[i].re - t.re;
					x[k].im = x[i].im - t.im;

					x[i].re += t.re;
					x[i].im += t.im;
				}
			}
		}

		for (int i = 1; i < Size - 1; i += 2) {
			x[i].re = -x[i].re;
			x[i].im = -x[i].im;
		}
	}

	/**
	 * @param DirInv
	 *            - direction of transform
	 * 
	 * 
	 * */
	/*
	 * public void fft2mpdht( String SrcModPic, String SrcPhPic, String DstPic,
	 * int DirInv, double Scale, double Expand ) { if (Expand<0) { throw new
	 * IllegalArgumentException("Expand < 0"); }
	 * 
	 * int filepic = 1;//open( SrcModPic, 1); int phasepic = 1;//open( SrcPhPic,
	 * 1);
	 * 
	 * 
	 * int modhead[] = new int [4];// _read (filepic, modhead, 8); int
	 * phasehead[] = new int [4] ; //_read (phasepic, phasehead, 8);
	 * 
	 * if (modhead[2] != 32 || modhead[2] != 32) { throw new
	 * IllegalArgumentException("Wrong type of file, must be 32 bit"); } int m0
	 * = modhead[0]; int n0 = modhead[1]; int m1,n1,m,n,dm,dn;
	 * 
	 * int modCommLen = modhead[3]; int phaseCommLen = phasehead[3];
	 * 
	 * 
	 * m1=(int)(m0*(1+0.01*Expand)); n1=(int)(n0*(1+0.01*Expand));
	 * 
	 * int k = 1;
	 * 
	 * while (k<m1) { k *= 2; }
	 * 
	 * m=k; dm=(m-m0)/2;
	 * 
	 * k = 2; while (k<n1) { k *= 2; }
	 * 
	 * n=k; dn=(n-n0)/2;
	 * 
	 * int reshead[] = new int [4];
	 * 
	 * reshead[0] = m; reshead[1] = n; reshead[2] = 64; reshead[3] = 0;
	 * 
	 * String comment = new String();
	 * 
	 * int buflen = (n>m) ? n : m;
	 * 
	 * Compl strin[]= new Compl[buflen]; Compl struse[]= new Compl[buflen];
	 * 
	 * int i,j,S=n/2;
	 * 
	 * Compl[] tab = new Compl[S+1];
	 * 
	 * 
	 * for(j=0;j<=S;j++) { tab[j].re=Math.cos(Math.PI*j/S);
	 * tab[j].im=-DirInv*Math.sin(Math.PI*j/S); }
	 * 
	 * 
	 * 
	 * float sfile[]=new float [n]; float ph[]=new float [n];
	 * 
	 * if (Scale==0) { Scale = (DirInv==1) ? 1 : 1.0/m/n; }
	 * 
	 * int desk=1;//open(DstPic,772,128);
	 * 
	 * long FreeMem = 1;//own_farcoreleft(), long pp = FreeMem/(8*n);
	 * 
	 * 
	 * int p,kluch=0;
	 * 
	 * if(pp>m0) { p = m0;} else { p = (int)pp; }
	 * 
	 * if(pp>=m) { kluch=1; }
	 * 
	 * int p_last,npr,p0;
	 * 
	 * Compl buf[]; Compl ptr[];
	 * 
	 * if(kluch==0) { p_last = m0%p; npr=(m0+p-1)/p;
	 * 
	 * buf =new Compl [n*8*p]; } else { buf =new Compl [n*8*m]; }
	 * 
	 * if(kluch==1){ BM (); return; //FIRST GOTO }
	 * 
	 * 
	 * 
	 * 
	 * for(i=0;i<dm;i++) { own_write( desk, strin, n*8); } p0=p;
	 * for(k=0;k<npr;k++){ if( k==npr-1 && p_last!=0 ){ p=p_last; }
	 * for(i=0,ptr=buf;i<p;i++,ptr+=n){ own_setmem(sfile,4*n,0);
	 * own_read(filepic, sfile+dn, n0*4); own_read(phasepic, ph+dn, n0*4);
	 * for(j=0;j<n;j++){ struse[j].re=Scale*sfile[j]*Math.cos(ph[j]);
	 * struse[j].im=Scale*sfile[j]*Math.sin(ph[j]); } fftdhts(struse, n, tab);
	 * for(j=0;j<n;j++){ ptr[j].re=struse[j].re; ptr[j].im=struse[j].im; } }
	 * own_write(desk,buf,n*8*p); }
	 * 
	 * own_setmem(strin,8*n,0);
	 * 
	 * for(i=m0+dm;i<m;i++){ own_write( desk, strin, n*8); } own_farfree(buf);
	 * COMM(); return; }
	 * 
	 * private void BM() { for(int i=0,ptr=buf;i<m;i++,ptr+=n){ if(i<dm ||
	 * i>=m0+dm){ own_setmem(ptr,8*n,0); continue; } own_setmem(sfile,4*n,0);
	 * own_read(filepic, sfile+dn, n0*4); own_read(phasepic, ph+dn, n0*4);
	 * for(int j=0;j<n;j++) { struse[j].re=Scale*sfile[j]*cos(ph[j]);
	 * struse[j].im=Scale*sfile[j]*sin(ph[j]); } fftdhts(struse, n, tab);
	 * for(int j=0;j<n;j++){ ptr[j].re=struse[j].re; ptr[j].im=struse[j].im; } }
	 * COMM(); }
	 * 
	 * private void COMM () {
	 * 
	 * own_farfree((void far *)sfile); own_farfree((void far *)ph);
	 * 
	 * _read( filepic, comment, modCommLen); _read( phasepic, comment +
	 * modCommLen, phaseCommLen );
	 * 
	 * close(filepic); close(phasepic);
	 * 
	 * comment[ modCommLen + phaseCommLen ] = 0; strcat( comment,
	 * "FFT2MP{(Module,Phase):");
	 * 
	 * if (DirInv==1) { strcat(comment,"Direct}"); } else{ strcat( comment,
	 * "Inverse}"); }
	 * 
	 * int cl = strlen(comment);
	 * 
	 * close(filepic); close(phasepic);
	 * 
	 * if(m==n) { EQ(); return; }
	 * 
	 * own_farfree(tab);
	 * 
	 * S=m/2; tab=new Compl [S+1];
	 * 
	 * for(int j=0;j<=S;j++){ tab[j].re=Math.cos(Math.PI*j/S);
	 * tab[j].im=-DirInv*Math.sin(Math.PI*j/S); } EQ(); }
	 * 
	 * private void EQ(){ long ind; if(kluch==1) { BM1(); return; } FreeMem =
	 * own_farcoreleft(); pp = FreeMem/(8*m);
	 * 
	 * if(pp>n) { p = n;} else { p = (int)pp; }
	 * 
	 * p0=p; p_last = n%p; npr=(n+p-1)/p;
	 * 
	 * buf = new Compl[m*8*p];
	 * 
	 * for (i=0;i<npr;i++) { Compl[] pointer = buf; lseek (desk, i*p0*8+8,
	 * SEEK_SET);
	 * 
	 * if( i==npr-1 && p_last!=0 ) { p=p_last; }
	 * 
	 * for( int j=0; j<m; j++, pointer +=p ) { own_read( desk, strin, 8*p);
	 * own_movedata(strin,pointer,8*p); lseek( desk, 8*(n-p), SEEK_CUR); }
	 * 
	 * for ( j=0; j<p; j++) { for ( k=0,ind=j; k<m; k++,ind+=p){
	 * strin[k]=buf[ind]; } for(k=0;k<m;k++){ struse[k].re=strin[k].re;
	 * struse[k].im=strin[k].im; }
	 * 
	 * fftdhts(struse,m,tab);
	 * 
	 * for(k=0;k<m;k++){ strin[k].re=struse[k].re; strin[k].im=struse[k].im; }
	 * for ( k=0,ind=j; k<m; k++,ind+=p) { buf[ind].re=strin[k].re;
	 * buf[ind].im=strin[k].im; } } lseek(desk, (long)i*p0*8+8, SEEK_SET);
	 * pointer=buf;
	 * 
	 * for(j=0; j<m; j++, pointer +=p ) { own_movedata(pointer,strin,8*p);
	 * own_write(desk, strin, 8*p); lseek(desk, 8*(n-p), SEEK_CUR); } } SVOB();
	 * }
	 * 
	 * private void BM1 () { for(int j=0;j<n;j++){
	 * for(i=0,ind=j;i<m;i++,ind+=n){
	 * 
	 * struse[i].re=buf[ind].re; struse[i].im=buf[ind].im; }
	 * 
	 * fftdhts(struse,m,tab);
	 * 
	 * for(i=0,ind=j;i<m;i++,ind+=n){
	 * 
	 * buf[ind].re=struse[i].re; buf[ind].im=struse[i].im; } }
	 * own_write(desk,buf,n*8*m); SVOB(); }
	 * 
	 * private void SVOB (){ own_farfree(buf); own_farfree(struse);
	 * own_farfree(tab); own_farfree(strin);
	 * 
	 * lseek(desk,6,SEEK_SET); _write(desk,&cl,2); lseek(desk,m*n*8,SEEK_CUR);
	 * _write(desk,comment,cl); close(desk); delete comment; return NoError; }
	 */
}
